Von-Neumann's and related scaling laws in Rock-Paper-Scissors type models 
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We introduce a family of Rock-Paper-Scissors type models with Z^ symmetry (A'' is the number 
of species) and we show that it has a very rich structure with many completely different phases. 
We study realizations which lead to the formation of domains, where individuals of one or more 
species coexist, separated by interfaces whose (average) dynamics is curvature driven. This type of 
behavior, which might be relevant for the development of biological complexity, leads to an interface 
network evolution and pattern formation similar to the ones of several other nonlinear systems in 
condensed matter and cosmology. 



The mechanisms leading to the enormous biodiversity 
observed in nature are still not fully understood. Rock- 
Paper-Scissors (RPS) type models incorporate some of 
the crucial ingredients associated with the dynamics of 
a network of competing species and they have been used 
as a powerful tool in the study of complex biological 
systems. In its simplest form, the RPS game describes 
the evolution of 3 species which cyclically dominate each 
other [iHSj (see also [H ^ for the pioneer work by Lotka 
and Volterra). The basic interactions are Motion, Re- 
production and Predation but generalizations, incorpo- 
rating more than 3 species and/or new interactions be- 
tween them, have been proposed in the literature [6til4j. 
Particularly, in |B] segregation processes and phase tran- 
sitions in predator-prey models with an even number of 
species have been investigated. 

RPS type models naturally lead to the formation of 
complex spatial patterns observed in some biological sys- 
tems p!]. Complex spatial structures also arise in many 
other systems. For example, interfaces in ideal soap 
froths and grain growth have a velocity v proportional 
to the mean curvature k at each point, which is at the 
core of the Von Neumann's law jTS] for the evolution of 
the area of individual domains. The evolution with time 
t of the characteristic scale L of such interface networks 
obeys the scaling law L ex i^/^, leading to the formation 
of cellular patterns of increasing size |T6H23] . Despite the 
different dynamical scaling laws, the evolution of inter- 
face networks in a cosmological context may also gener- 
ate similar spatial patterns to the ones observed in soap 
froths [SUES]- 

Consider a family of models where individuals of vari- 
ous species are distributed on a square lattice of size M 
at some initial time. The different species are labelled 
by the number i (or j) with i,j — 1, ..., A^, and we make 
the cyclic identification i — i -\- k N where k is an inte- 
ger. The number of individuals of the species i will be 
denoted by It. In addition to individuals, there are also 
empty spaces which shall be denoted by E. In this pa- 



per, except if stated otherwise, we shall assume that the 
initial distribution is random and that the number den- 
sities of the various species rii = li/M are all identical 
at the initial time. The number density of empty spaces 
is given by tle = IeI-N' and its initial value is equal to 
0.1 in all the network simulations described in this paper. 
At each time step a random individual (active) interacts 
with one of its four nearest neighbors (passive) . The unit 
of time At = 1 is defined as the time necessary for M in- 
teractions to occur (one generation time). The possible 
interactions are classified as Motion 



Reproduction 



iE -^ Ei , or ij — )■ ji . 



iE ^f ii . 



or Predation 

i{i ^ a) ^ lE , or j(i -\- a) ^ iE , 

where a = 1, ...,amax with amax — N/2 if A^ is even 
or Umax = (A^ — l)/2 if A^ is odd. We shall denote the 
corresponding probabilities by m^ (Motion), r^ (Repro- 
duction), pna (left-handed Predation) and pma (right- 
handed Predation) . This family of models has a Z^ sym- 
metry if m, = m,r, = r, p^i^ = PLa and pma = PRa for 
all i. Figure 1 represents the different predation inter- 
actions in a model with 5 species having Z^ symmetry. 
Throughout this paper, we shall assume that Af — 600^ 
and that the Zpf symmetry is preserved with the follow- 
ing interacting probabilities: to = 0.5, r = 0.25 (or zero, 
if the passive is not E), p = 0.25 (or zero, if the passive 
is E or if the passive is not a prey of the active). The 
Zn symmetry is, in general, not associated with curva- 
ture driven dynamics. For example, the standard RPS 
model has a Z^ symmetry but the dynamics of the spa- 
tial patterns is not curvature driven in this model. In 
fact, we shall show that the dynamics is curvature driven 
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FIG. 2: (Color online) Snapshot of the evolution of model 
I. The 2 species are represented using different colors on the 
left panel. On the right panel the black dots represent the 
distribution of empty spaces. 



FIG. 1: (Color online) Predation interactions in a model with 
5 species having Z^ symmetry. 



only in realizations which result in the formation of do- 
mains, where individuals of one or more species coexist, 
separated by interfaces whose dynamics is controlled by 
interactions of equal strength between competing species. 

Let us start by introducing a simple model with N = 2 
having p^i = pri ^ (model I). This model has 2 com- 
peting species which tend to distribute themselves into 
separate domains bounded by interfaces where most of 
the action occurs. Predation happens mainly at the in- 
terfaces whose thickness is directly related to the mobil- 
ity. This is clearly shown on Fig. 2 where a snapshot of 
the network evolution of model I is presented. On the 
left panel each species is represented by a different color 
while on the right panel the black dots represent the dis- 
tribution of empty spaces. Fig. 2 clearly shows that the 
empty spaces, which are a result of Predation between 
individuals of the competing species, are located at the 
domain borders. The average thickness e of the interfaces 
does not change with time. 

For an interface with curvature radius p and thick- 
ness e ^ p, the average number of attacks per unit time 
from individuals outside the border is proportional to 
the outer interface length (proportional to p + e/2) while 
the average number of attacks from individuals inside 
the border is proportional to the inner interface length 
(proportional to p — e/2). The difference between the 
average number of attacks per unit of time from outside 
and inside the border is proportional to e. On the other 
hand, the number of attacks necessary to modify the do- 
main radius by Ap ^ p is proportional to the interface 
length (ex p). This implies that the value of the velocity 
of the interface is on average proportional to its curva- 
ture {v = Ck — C p^^ , where C is a positive constant), 
which is typical of non-relativistic interfaces in condensed 
matter ,25j . 

The average evolution of the area of a compact simply 
connected domain with no vertices is then given by d = 
§ vdl = —C § ndl = — 27rC, where dl is the infinitesimal 



interface arc length and a dot denotes a derivative with 
respect to t. Note that the domain area can be calculated 
at any given time by counting the number of individuals 
inside the domain. If the domain is compact, but not 
necessarily simply connected, then 



\9] 



27rC(g-l), 



(1) 



where argi represents the area of a compact domain with 
genus g and no vertices. Eq. (fTl) implies that the area de- 
creases (increases) with time depending on whether g = Q 
{g > 0). The genus dependency accounts for the contri- 
bution of the decrease of the area of each hole to the 
growth of the area of the domain, li g = then, on 
average, the evolution of the area with time is given by 



a{t) = a(0) (1 - t/Q , 



(2) 



where tj, is the time of collapse. The average time evo- 
lution of the area a(t) of an initial circular domain with 
f; = is illustrated in Fig. 3. The domain area a(t) was 
calculated by counting the number of individuals of the 
species inside the domain. We verified that a nearly iden- 
tical result is obtained using the relation a oa I^. The 
solid red line represents the average result in an ensemble 
of 43 simulations using model I and the dashed blue line 
shows the theoretical evolution given by Eq. (pi), with 
tc calculated as the median collapse time. Fig. 3 shows 
that the agreement between the analytical and numerical 
results is very good. 

In the case of an interface network without junctions 
the evolution of the total number of domains, Nd, with 
time can be obtained using Eq. (flj). If the fraction /[qj 
of the total number of domains with genus (7 = is a 
constant (/[qj = 0) then [TB] 



Nd 



-f[o] 



(3) 



Here digi ~ ^[g]d is the average area of a compact domain 
with genus g, a — A/Nd is the average domain area, 
A = X]n=o ^[g] ^^ ^^® area of the entire system and A^g] 
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FIG. 3: (Color online) The time evolution of the area of a 
circular domain. The solid red line represents the average 
result in an ensemble of 43 simulations using model I and the 
dashed blue line is the theoretical prediction. 



is the total area occupied by domains with genus g 
both /[Qi and htQ^ are time independent, then Njj oc 
The characteristic length of the network defined by L 
{A/NdY^"^ evolves as 



If 



Nl. 



L oc N^^'^ ex <i/2 , 



(4) 



Fig. 4 shows the evolution of the number of empty 
spaces Ie with time t for model I. The cyan dots (light 
grey in black and white) represent the results of 20 inter- 
face network simulations while the solid line represents 
the average. The scaling exponent A defined by Ie oc t^ 
is A = —0.450 ± 0.027 at one-sigma. The number of do- 
mains No is equal to the ratio between the total area 
A and the average domain area L^. It is also propor- 
tional to the ratio between the total interface length Lt 
and the average domain perimeter which is proportional 
to L. Hence, A/L'^ ex Lt/L implying that L ex A/Lt- 
On the other hand, the average thickness e of the inter- 
faces does not change with time and consequently the 
total interface length L^ is proportional to Ie- Hence, L 
is inversely proportional to the number of empty spaces 
{L (X I]^ ) and consequently one would expect an av- 
erage scaling solution with Ie oc t^, with A = —0.5. 
The fact that the numerical value of A is very close to 
the theoretical one demonstrates that the interface net- 
work evolution is already attaining the expected scaling 
regime. The increase with time of the dispersion between 
the values of Ie for the various simulations is associated 
with the growth of the characteristic length scale of the 
network (see [26] for a detailed analytical discussion). 

li V = Ck then the evolution of the area of a single 
domain with genus 5 = and £ vertices is given by 



hi = — C (p ndl = — C 




(5) 



where Op represent each of the i discontinuous angle 
changes at the vertices. If the interface network has only 



FIG. 4: (Color online) Evolution of the number of empty 
spaces Ie with time t for model I. The cyan dots (light grey in 
black and white) represent the results of 20 interface network 
simulations while the solid line represents the average. The 
scaling exponent A defined by Ie "' *^ 
-0.450 ±0.027. 



oc i is estimated as A 



Y-type junctions which meet at an angle of 27r/3, then 
one obtains the Von Neumann's law [TS] 



ae = -C 



27r - £ TT - 



27r 



= C-{e-6), (6) 



implying that the area domains with i < 6 (i > 6) de- 
creases (increases) proportionally to time. The evolution 
of interface networks in RPS type models is not deter- 
ministic and consequently the Von Neumann's law can 
only be valid on average. Fig. 5 shows the evolution of a 
Y-type junction in model IV where pLa = PRa for all a 
(tliis model will be described later in more detail), start- 
ing from the initial configuration shown in the left panel. 
The right panel represents the most frequent species at 
each lattice point from i = 1 x 10^ to t = 2.5 x 10^, show- 
ing that the average angles at the vertex all tend to an 
average value of 27r/3. 

The evolution of the total number of domains Ne with 
time, in the case of an interface network with Y-type 
junctions, can be obtained using Eq. ^. If the frac- 
tion fe of the total number of domains with i edges is a 
constant {fg = 0) then 



Ne 

Nr 



> —We 

^ — ^ \ n /I 



(7) 



Here ai = hio, is the average area of a domain with 
£ edges, a = A/Ne is the average domain area, A = 
'^i=i ^f i*^ the area of the entire system, Ai is the total 
area occupied by domains with £ edges and the function 
wi accounts for the fact that the collapse of an individual 
domain might lead to the merger of some of the surround- 
ing domains (wg = 1 if all domains are of different types 
and wg > 1 otherwise). If the interface network is in a 
scaling regime with time- independent /^, he and we, then 
again Ne oc —Nf^. Consequently, the scaling law given 




FIG. 5: (Color online) Evolution of a Y-type junction in 
model IV where pLa ~ Prc for all a. The left panel shows the 
initial configuration with different species represented by dif- 
ferent colours. The right panel represents the most frequent 
species at each point from t = 1 x 10^ to t = 2.5 x 10^ (note 
that the average angles at the vertex are all approximately 
equal to 27r/3). 




FIG. 6: (Color online) Snapshots of 4 different 600^ simula- 
tions with N = 6 after 4000 generations. The different panels 
show the results obtained using model II (upper left), model 
III (upper right), model IV (lower left) and model V (lower 
right). 



in Eq. ^ also applies to interface networks with Y-type 
junctions. 

Now we shall demonstrate that the curvature driven 
interface dynamics of model I is common in RPS type 
models. Let us consider a model with iV = 6. By tak- 
ing Pli = P and pm = pl2 = Pr2 = Pls = PR3 = 
(model II) one ensures that the domains which appear 
in the simulations are populated with one of the follow- 
ing non-interacting species partnerships: {1,3,5} and 
{2,4,6}. The interactions (of equal strength) between 
the 2 different partnership domains occur mainly at the 



border where a large number of empty spaces are contin- 
uously being generated (see video [27| and Fig. 6 (upper 
left)). If one considers the case with p^i — p]^2 — Pi 
Pri — PR2 = PL3 = PR3 = (model III) then there are 3 
possible partnership domains containing non-interacting 
species: {1, 4}, {2, 5} and {3, 6} thus leading to an inter- 
face network with Y-type junctions (see video [55] and 
Fig. 6 (upper right)). In the case with p^a = PRa = P 
for a = 1,2,3 (model IV), there are no partnerships 
since every species is linked in a bidirectional way to all 
other species thus generating an interface between any 2 
given domains. This gives rise to an interface network 
with Y-type junctions and 6 different domain types (see 
video [29] and Fig. 6 (lower left)). If one now takes 
Pli = Pri = Pl2 = P and pb.2 = Pl3 = Pr3 = (model 
V) then there are 2 possible domains defined by {1, 3, 5} 
and {2, 4, 6} but due to the non-zero unidirectional prob- 
ability pr2 of interaction between species in the same 
domain, spiral patterns do form (see video |30) and Fig. 
6 (lower right)). Inside the 2 different domains the in- 
teractions are those of the standard RPS model with 3 
species. In the boundary there is an (average) equilib- 
rium between the predation probabilities from individu- 
als on either side of the wall (for example, individuals of 
the species 1 predate and are predated by individuals of 
the species 2 and 6) . In models IV and V the colors light 
blue, dark blue, red, maroon, green and yellow represent 
species 1 to 6, respectively. This study represents a sig- 
nificant extension with respect to previous work and, to 
the best of our knowledge, it is the first time that models 
III, IV and V have been studied in the literature. Note 
that model V is very different from the 6 species exten- 
sion of the 5 species Rock-Paper-Scissors-Lizard-Spock 
model proposed in [T^], defined by pm = pl2 — P and 
Pli = PR2 = PR3 = PL3 = (see video [SI]), whose 
dynamics is not curvature driven. 

The scaling parameters A calculated numerically for 
aU models are: A = -0.450 ± 0.027 (model I), A = 
-0.467±0.034 (model II), A = -0.465±0.028 (model III), 
A = -0.421 ± 0.018 (model IV) and A = -0.429 ± 0.029 
(model V). In the case of model V, the empty spaces 
also appear associated to the spiral patterns. Hence, the 
parameter A for this model was obtained using only the 
empty spaces which have as some of the four immediate 
neighbors individuals from the 2 groups: {1,3,5} and 
{2,4,6}. In all cases the scaling parameter A is reason- 
ably close to the expected value A = —0.5. The devia- 
tions can be attributed to the finite size and dynamical 
range of the simulations. We have verified that the re- 
sults obtained for models I, II, III, and IV are weakly 
dependent on the value of m, as long as the thickness 
of the interfaces is much smaller than the box size. In 
the case of model V we have found that the network 
evolution is curvature dominated only if ?7i > 0.5. For 
smaller values of m there are other effects which have a 
significant impact on the network dynamics, such as lo- 
cal partnerships along the borders of the domains [5] (for 
example, between species 1 and 4 which do not predate 



each other). These effects are outside the scope of the 
present paper and will studied in detail in future work. 
In summary, we have shown that curvature driven in- 
terface dynamics, analogous to that observed in other 
physical systems in condensed matter and cosmology, is 
common in the family of RPS type models investigated 



in this paper and may be crucial to the understanding of 
biological complexity. 
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